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INTERPOLATION ERROR ESTIMATES FOR HARMONIC COORDINATES ON 

POLYTOPES 

ANDREW GILLETTE* AND ALEXANDER RANDt 


Abstract. Interpolation error estimates in terms of geometric quality measures are established for harmonic coordinates 
on polytopes in two and three dimensions. First we derive interpolation error estimates over convex polygons that depend 
on the geometric quality of the triangles in the constrained Delaunay triangulation of the polygon. This characterization is 
sharp in the sense that families of polygons with poor quality triangles in their constrained Delaunay triangulations are shown 
to produce large error when interpolating a basic quadratic function. Non-convex polygons exhibit a similar limitation; large 
constrained Delaunay triangles caused by vertices approaching a non-adjacent edge also lead to large interpolation error. While 
this relationship is generalized to convex polyhedra in three dimensions, the possibility of sliver tetrahedra in the constrained 
Delaunay triangulation prevent the analogous estimate from sharply reflecting the actual interpolation error. Non-convex 
polyhedra are shown to be fundamentally different through an example of a family of polyhedra containing vertices which are 
arbitrarily close to non-adjacent faces yet the interpolation error remains bounded. 


1. Introduction. Interpolation error estimates over triangles or higher-dimensional simplices are an 
essential component of the most common finite element analyses. The simplest approach involves restrict¬ 
ing attention to simplices of bounded aspect ratio, which, in two-dimensions, is equivalent to the familiar 
minimum angle condition. However, this is an overly restrictive approach and identical estimates have been 
established over broader classes of simplices. In two-dimensions, improved estimates involve the maximum 
angle condition [7, 51, 42] which can be generalized to higher dimensions [45, 52, 76]. These types of condi¬ 
tions and ideas can then be applied to estimates in different norms [83, 1, 70], other types of elements [31, 2] 
or, in some sense, even to finite volume methods [65]. The interpolation estimates are sharp in the sense that 
there exist sequences of function/simplex pairs that realize the estimate [7, 5]. While the maximum angle 
condition is not strictly necessary for convergence of finite element methods [43] , the known counterexamples 
involve finite element spaces which contain subspaces corresponding to shape-regular meshes. Thus, some 
form of shape regularity still appears to be a necessary ingredient for successful interpolation. 

Polygonal finite elements can be established as a generalization of triangular finite elements by attempting 
to preserve many properties from the triangular case. This essential framework has been applied to a variety 
of different finite element contexts [81, 35, 84, 91, 64] and shares much of the underlying theory with related 
numerical methods for polygonal meshes: mimetic methods [4, 18, 61, 19, 12], virtual element methods [8, 
10, 9, 21, 62], weak Galerkin methods [88, 87, 66], compatible discretization operator schemes [13, 14, 55], 
and some “meshfree” methods [58, 82, 6]. 

Generalized barycentric coordinates (GBGs) are a common approach to defining basis functions for 
polygonal and polyhedral finite element methods. A set of GBGs on an n-dimensional polytope will include 
in their span the set of linear functions on K", but the GBGs themselves are in general not polynomials. 
Gonstruction of GBGs is not unique and many competing variants have been developed. These include 
Wachspress coordinates [86], mean value coordinates [38], Sibson coordinates [79], maximum entropy coor¬ 
dinates [80], Poisson coordinates [60], and harmonic coordinates [46, 64]. Motivation for GBGs include data 
interpolation [79], computer graphics [46], and mesh generation [20], as well as finite element applications. 
Interpolation error, the focus of the paper, requires the development uniform estimates over some class of 
polytopes for a particular generalized barycentric construction. Since all GBGs do not admit estimates 
over the same geometric restrictions, identification of coordinates which admit the widest set of polytopes 
provides one method of differentiating particular constructions. 

Among all possible GBGs, harmonic coordinates are optimal in the context of interpolation error: in¬ 
terpolation by harmonic coordinates minimizes the i7^-seminorm over all functions satisfying the requisite 
boundary data. Thus, a bound on the interpolation error of harmonic coordinates in terms of geometric 
properties of the domain serves as a benchmark for the error analysis of all other types of GBGs. Generalizing 
geometric conditions from simplices to polytopes, however, must be done with some care; some conditions 
which are equivalent on triangles are no longer equivalent on polygons. For example, requiring a minimum 
angle is equivalent to bounded aspect ratio on triangles, while a generic convex polygon can have no small 
angles but poor aspect ratio, as in the case of a rectangle with a large ratio of width to height. Prior error 
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analysis is primarily aimed at generalizing the minimum angle condition on triangles. For harmonic coor¬ 
dinates error estimates were established for convex polygons with bounded aspect ratio [41], and for other 
GBCs, some additional restrictions are were placed [41, 71]. In dimension larger than two, similar estimates 
have been established for Wachspress coordinates for some families of shape-regular polytopes [39]. 



Fig. 1.1: Three types of triangles from the perspective of interpolation quality metrics: (A) a good aspect 
ratio triangle, (B) a poor aspect ratio triangle with a single small angle (and thus modest circumradius) and 
(C) a poor aspect ratio triangle with a large angle (and thus large circumradius). Standard interpolation 
error estimates hold for the first two types but fail for the third. 


In this paper, we broaden the analysis of harmonic coordinates from [41] to non-bounded aspect ratio 
polygons and to three-dimensional polyhedra. We connect error in the harmonic coordinates to errors in a 
piecewise linear approximation for certain triangulations/tetrahedralizations of the polytope. Despite the 
fact that piecewise linear interpolation with respect to a triangulation of the domain and interpolation by 
harmonic GBCs appear to be at opposite ends of the accuracy spectrum in our previous error analysis, we 
show here that for convex polygons, the constrained Delaunay triangulation provides an interpolation error 
that is essentially as good as the harmonic coordinates. For non-convex polygons, we establish a similar re¬ 
sult connecting large circumradius triangles in the constrained Delaunay triangulation to large interpolation 
errors, although this requires few additional assumptions to avoid some other degenerate configurations that 
are only possible in the non-convex setting. We also analyze harmonic coordinates in three dimensions and 
establish error estimates for bounded aspect ratio polyhedra. Features of polyhedra that impact interpola¬ 
tion quality are described in both the convex and non-convex cases. Interpolation on convex polyhedra is 
analogous to convex polygons although due to sliver tetrahedra, the connection to the Delaunay tetrahedral- 
ization cannot be made. Rather, poor quality is characterized by vertices being near non-adjacent faces. For 
non-convex polyhedra, this characterization does not hold and we give an example demonstrating degenerate 
limiting geometry with bounded interpolation error. 

Section 2 begins by stating sharp interpolation results for triangles as a fundamental building block 
of the analysis. Section 3 introduces the key ideas and definitions of generalized barycentric interpolation. 
Sections 4 and 5 contain analysis of harmonic coordinates in two and three dimensions, respectively. Ap¬ 
pendices A, B, G, D, and E contain a few well-known results used in our analysis from the theory of Sobolev 
spaces, partial differential equations, computational geometry, and numerical analysis. 

2. Barycentric Interpolation on Triangles. The classical a priori analysis of the finite element 
method involves two steps. First the finite element error is bounded (up to some constant factor) by the 
error in the best possible interpolant in a particular finite element space; this result is often called Cea’s 
Lemma. Second, a specific interpolant is constructed and is demonstrated to have a suitably small error in 
terms of some geometric quantities related to the mesh. 

The most common interpolant selected for linear triangular elements is derived from barycentric coor¬ 
dinates. The barycentric functions for a triangle are the three affine functions taking value 1 at one 

vertex of the triangle and taking value 0 at the other two vertices. The standard Lagrange interpolant on 
triangle T is defined in terms of the barycentric coordinates by 
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Fig. 1.2: A number of convex polygons with different shape quality characteristics, some of which are distinct 
from the three classes of triangles: (A) a polygon with modest aspect ratio and no short edges; (B) a very 
narrow polygon that does not contain any short edges or large angles in the Delaunay triangulation; (C) a 
polygon with a large angle in the constrained Delaunay triangulation; (D) a polygon with good aspect ratio 
but some very short edges and small angles in the Delaunay triangulation. The Delaunay triangulation of 
each polygon is shown. (Note: since the polygons are convex, no boundary edges need to be constrained.) 
Even if it is part of a larger linear segment (as occurs in (B), (C), and (D)), each segment between vertices 
shown is treated as an independent edge for the purpose of the upcoming construction (3.2). 


for any function it : T —> M that is smooth enough to admit pointwise values. In the construction of finite 
element spaces, this yields two key properties. First, the interpolant along a boundary edge of a triangle only 
depends on the values at the two vertices of the triangle, which allows these functions to be stitched together 
continuously over the entire mesh. Second, the space spanned by the barycentric coordinates contains all 
linear polynomials, which is a crucial piece of the error analysis. 

The simplest analysis of the finite element method (widely adopted in the initial exposition of the 
method in many texts [25, 17, 34, 92]) involves interpolation estimates on triangles under the minimum 
angle condition: for minimum angle 0 > 0, there exists a constant Cq such that for all triangles T with 
minimum angle larger than Q and all functions u G H^{T), 

(2.2) jit ^ Cghx . 

The quantity hx indicates the length of the longest edge of a triangle T. Under the minimum angle condition, 
this quantity is equivalent (up to a constant factor) to many other measures of the triangle’s “size,” such as 
shortest edge length, inscribed radius, circumradius, diameter, etc. 

Nevertheless, some triangles with arbitrarily small minimum angles exhibit small interpolation errors. 
Put differently, only triangles with large angles (near 180°) yield large interpolation errors. Figure 1.1 
depicts the three different classes needed in the rehned analysis. Under the maximum angle condition, the a 
priori estimate (2.2) can be established [83, 7, 45], again using the longest edge of the triangle to quantify 
a triangle’s size, which is no longer equivalent to all of the other measures. 

Further, a uniform estimate can be established for all triangles that captures the dependence of the 
constant on the largest angle. The most natural form of this estimate avoids any explicit reference to 
triangle angles by quantifying triangle size using the circumradius. 

Theorem 2.1. There exists a constant C2.3 such that for all triangles T and all functions u G H^(T), 

(2.3) |ii — — C 2 . 3 RT \'^\h‘^{t) > 

where Rt denotes the circumradius of triangle T. Throughout the paper, we will use the notation Rt 
to denote the circumradius of triangle T; for an example, see Figure D.2 in the appendix. While some 
authors have directly connected error estimates over triangles to the circumradius [52, 69, 16, 49, 50], this 
estimate is also stated in several other essentially equivalent ways in the literature, such as hxj cos(q:t/ 2 ) 
or hx/ s\n{aT), where a is the largest angle in triangle T [45, 42]. 
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3. Interpolation with Generalized Barycentric Coordinates. Generalized barycentric coordi¬ 
nates extend some key properties of barycentric coordinates to polygons in a way that allows for analogous 
finite element spaces to be constructed for polygonal meshes. Letting P be a polygon with vertices 
the functions are called generalized barycentric coordinates if they satisfy: 


(GBCl) 

(GBG2) 

(GBC3) 

(GBC4) 

(GBG5) 

(GBC6) 


Non-negativity: Ai > 0 on P; 

n 

Linear completeness: For any affine function L : P —>■ K, L = L{vi)Xi; 

i—1 

Invariance: If P : — >■ is the composition of rotation, translation and/or uniform scaling 

operations, then Ai(x) = Af (P(x)), where A^ denotes the barycentric coordinate on P; 

n 

Partition of unity: = 1; 

n 

Linear precision: ^ViAi(x) = x; 

i=l 

Interpolation: Xi { vj ) = Sij . 


These properties are not independent axioms. In particular, the last three properties can be derived from 
the first two. We state them together as they identify the complete set of standard properties of generalized 
barycentric coordinates. 

While standard barycentric coordinates on a triangle are unique affine functions that depend only on the 
location of the vertices, the properties (GBG1)-(GBG6) do not define a unique set of generalized barycentric 
coordinates in general. This non-uniqueness has allowed the development of various generalized barycentric 
coordinate functions tailored to distinct application contexts. Gomparison of the different coordinate types 
centers around three distinguishing factors: constraints on the polygonal domain P, ease of computation of 
the coordinates, and the smoothness of the resulting interpolant. The ‘interpolant’ associated to a set of 
generalized barycentric coordinates is analogous to the definition in (2.1), with the sum being taken over 
all the vertices in P. From this standpoint, we can identify a spectrum of coordinate types ranging from 
simplicity of implementation to simplicity of analysis. On one end lie “triangulation” coordinates and on 
the other lie harmonic coordinates, while other types lie somewhere in between. 

Triangulation coordinates are the most straightforward generalized barycentric coordinates both con¬ 
ceptually and from the standpoint of computation. In this case, the polygon is triangulated and a piecewise 
function is constructed with the barycentric coordinates over the triangles. This approach comes with some 
notable downsides. The triangulation coordinates are not smooth on the interior of the polygon: the inter¬ 
polant is in H^{P) but has discontinuous derivatives along the edges of the triangulation. Moreover, a poor 
quality triangle can cause large interpolation errors, even if the polygon satisfies certain quality metrics that 
would make it suitable for other generalized barycentric constructions. Hence, triangulation coordinates are 
sensitive to the triangulation selected. 

Harmonic coordinates occupy the opposite end of the spectrum: they are computationally expensive 
(as they in general have no explicit formula) but are optimal from the perspective of smoothness of the 
interpolant. Specifically, the harmonic coordinates are constructed as the solution of the partial differential 
equation: 


f AXi = 0, on P, 

\ X^ = g^. on dP. 


The boundary condition gi : dVl —>■ K is the piecewise linear function satisfying 


(3.2) = A) linear on each edge of H. 

Since this is a linear PDF, the resulting interpolant can also be characterized as a solution to the differential 
equation: 


(3.3) 


A {Ipu) 

IpU 
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0 on P, 
gu on dP. 


Here, is the piecewise linear function equal to u at the vertices of P. We use the same notation for harmonic 
coordinates as the barycentric interpolant on triangles without confusion, since harmonic coordinates on 
triangles produce the standard linear interpolant. 

These coordinates are optimal in the sense that they minimize the norm of the gradient of the interpolant 
over all functions satisfying the boundary conditions, that is, 

(3.4) Ipu = argmin | : v = guOn 9P|. 

This property, called Dirichlet’s Principle (see Appendix B), makes harmonic coordinates ideal for finite 
element error analysis, since the key estimate in the typical analysis relies on showing that the interpolation 
operator is bounded from into see [41] for more details. 

Many of the other generalized barycentric coordinate types have been developed with essentially the 
same goal of approximating the harmonic coordinates to avoid large gradients in the interpolant using a 
construction that does not require the solution of a PDE. While we only explicitly address harmonic coor¬ 
dinates in this paper, the characterization from (3.4) implies that any situation where harmonic coordinates 
produce large gradients applies to all generalized barycentric coordinates. 

4. Harmonic Coordinates on Polygons. 

4.1. Estimate Based on Triangulation Coordinates. We begin with an explicit bound on the 
harmonic coordinates’ interpolation error based on triangulation coordinates. This essentially formalizes an 
expected fact in light of (3.4): interpolation using harmonic coordinates is no worse than triangulating and 
using a piecewise linear interpolant on the triangulation. 

Theorem 4.1. There exists a constant (74.1 > 0 such that for any polygon P, possibly non-convex, all 
functions u G PI^{P) and all triangulations T of P, 

(4.1) \u- Ipu\pi(p^ < (74.1 ^max \u\p 2 ^p^ , 
where Rp is the circumradius of a triangle T. 

Proof. Without loss of generality, we restrict our analysis to a generic diameter one polygon P. Let Pu 
denote the linear approximate of u defined by the Bramble-Hilbert lemma; see Theorem C.l in Appendix C. 
Also noting that linear precision of the harmonic interpolant ensures that = IpPu, we begin by estimating 
the error with the triangle inequality: 

2 2 2 2 2 

(4.2) \u — Ipu\pi^p^ 'fi\u — Pu\p[i(^p'j + \Pu ~ ^P'^\h^{P) = Ir ~ P'u\h^(P) \^P (Pt* ~ ^)Ihi(P) ■ 

Next we appeal to Dirichlet’s principle (see Theorem B.l from Appendix B) and observe that the H^- 
seminorm of the harmonic interpolant is dominated by that of any piecewise linear interpolant that satisfies 
the same boundary conditions. Thus, for any triangulation T of P, 

(4.3) \Ip {pu — m)I^i(p) < |/r (Pu — '*^)Ihi(t) ■ 

T&r 


Combining (4.2) and (4.3) gives 

\u — Ipu\pi^p^ < |u — P«l^ri(p) + {Pu — '*^)lpi(T) ■ 

Ter 

Now while has been selected to be near u, it is not directly connected to Ipu. So we insert u into the 
estimate of that term and again apply the triangle inequality: 

|m — 7pu|pi(-p^ < I'M — Pii|pi(p) + ~ '^\m{T) + 1'*^ ~ ^Tu\pll(J.^^^ 

Ter 

= 2 ju — P-ulpl(p) + ^ ^ ~ ^T'U,\pif^p'^ ■ 

Ter 
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The first term is bounded by its construction from the Bramble-Hilbert lemma; note that the constant Cc.i 
below comes from Theorem C.l. Ij-u is the piecewise linear interpolant of u on triangle T and thus its 
interpolation error in the second term is bounded by Theorem 2.1. 


\u 




tgt 


< 


2 C'c 1 + C 2 3 max Rp 
tgT 




Since maxT^r Rt > 1/2 for any unit diameter polygon, (4.1) holds with C 4 ,i := iCc.i + C' 2 . 3 . □ 

The result just proved begs the question of how to find a triangulation that gives the sharpest estimate. 
For some polygons with interior angles near 180°, certain triangulations may yield arbitrarily poor interpo¬ 
lation properties even though other triangulations yield much smaller errors. See Figure 4.1 for an example. 
In attempt to find the best triangulation for a given polygon, a natural starting point is the Delaunay tri¬ 
angulation, which is based on a geometric criterion (the empty circumball property) that tends to avoid 
large circumradii triangles. Since we expect the boundary of the potentially non-convex polygon to be in the 
triangulation, we are required to use a constrained Delaunay triangulation [24, 77] , which is a generalization 
of the Delaunay triangulation that allows for required segments in the triangulation. 

Delaunay and constrained Delaunay triangulations have many optimality properties. Most well known, 
the Delaunay triangulation maximizes minimum angle [54, 78], maximizes mean inradius [53], and minimizes 
L'P error for the function jxj [27, 73, 22]. The most closely connected property is shown in [72, 67]: 
the Delaunay triangulation minimizes the H^-semi norm of the interpolant, but not necessarily the error. 
Further, [74] contains some explicit connections between Delaunay triangulation and the interpolation of 
harmonic functions. 



Fig. 4.1: Triangulations of a polygon that leads to a po 
The constrained Delaunay triangulation on the right i 



r (left) and good (right) interpolation error estimate, 
associated with bounded interpolation error. 


Finally, we emphasize that Theorem 4.1 applies to all planar polygons, including non-convex and non- 
bounded aspect ratio polygons that were not included in our original analysis of harmonic coordinates [41] . 
This is particularly useful as there are various kinds of polygon degeneracy that do not occur when only 
considering triangles. See Figure 1.2 for some examples. As a result, fully understanding interpolation 
estimates for harmonic coordinates on all polygons is reduced to two issues: ( 1 ) showing that the constrained 
Delaunay triangulation (CDT) for a particular polygon has maximum circumradius bounded by the polygon 
diameter, or (2) demonstrating that the error is in fact large when the CDT contains a large circumradius 
triangle. 

For convex polygons, we will show that this is an essentially complete way of characterizing the interpo¬ 
lation error of harmonic coordinates via the interpolation error of triangulation-based coordinates, excluding 
some specific families of degenerate polygons. Non-convex polygons lead to a similar conclusion, although 
the set of excluded degenerate families is somewhat larger in that setting. 

4.2. Sharpness: Convex Polygons. In certain situations, harmonic coordinates cannot improve upon 
interpolation via triangulation coordinates on the constrained Delaunay triangulation. To make this connec¬ 
tion explicit, we begin with a fact connecting large circumradii in the constrained Delaunay triangulation 
with obtuse triangles involving an edge along the boundary of the polygon. 

Lemma 4.2. Let P be a polygon and let i?* = vnaxT^T Rt, where T is the set of triangles in the 
constrained Delaunay triangulation of P. If R^. > diam(P), then there is a triangle & T such that (i) 
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Rt, = R* and (ii) the longest edge ofT^ is on the boundary of P. 

Proof. Suppose i?* > diam(P). Let T G T such that Rt = R*. If T is acute, then the circumcenter of 
T lies inside the triangle and thus diam(P) > i?*, a contradiction. So T is obtuse. Now suppose that the 
longest edge of T is not on the boundary of P and thus is not a constrained edge. In this case, the longest 
edge of T is shared by an adjacent triangle T' whose circumradius is at least as large as i?* by the Delaunay 
property; see Figure D.2 and Lemma D.l in Appendix D. Since T has a maximal circumradius by definition, 
the two circumradii must be the same. 

We can continue traversing the triangulation by this logic until we reach a triangle along the boundary 
of P. This walk will never return to a previously reached triangle, since (1) all triangles encountered must be 
obtuse and (2) two adjacent obtuse triangles cannot share the same longest edge in a Delaunay triangulation. 
Thus the length of the longest edge of the triangles in the walk is strictly increasing. Since there are a finite 
number of triangles, the walk must terminate in a triangle satisfying the conditions specified. □ 

Next, we demonstrate that the harmonic coordinate interpolant Ipu will not allow the standard inter¬ 
polation error estimate over any family of convex polygons with degeneracies that are essentially similar to 
a family of triangles with arbitrarily large angles. To this end, we consider a sequence of convex polygons, 
where one vertex approaches the interior of a non-incident edge and demonstrate that no constant can sat¬ 
isfy the error estimate. Let be a sequence of tuples where Pi is a convex polygons, is a 

vertex of Pi, and is an edge of Pi not incident to v^. We assume that this sequence satisfies the following 
assumptions: 

(AI) diam(Pi) = 1; 

(A2) Oi lies on the x-axis; 

(A3) Vi lies on the positive j/-axis; 

(A4) Pi lies below a line through with non-negative slope; 

(A5) There exists a c„ > 0 such that dist(o,Ve) > Cy, where is either endpoint of Ci and o is the origin; 

(A6) dist(vi,ei) —>■ 0 as i —>■ oo. 

Assumptions (AI)-(A4) can be taken without loss of generality via translation, scaling, rotation and 
reflection of the domain as needed. Note that after (AI)-(A3) have been accommodated, the convexity of Pi 
ensures that Pi lies below some line through Vi- If that line has negative slope, reflect Pi across the y-axis 
to accommodate assumption (A4). 

Assumptions (A5) and (A6) are the properties that ensure the polygons are degenerating in a way 
incompatible with the desired interpolation error estimates. For instance, the assumptions allow the case of 
a sequence of triangles with one angle approaching 180° (i.e. an unbounded circumradius) by setting Pi to 
be the triangle with vertices (—1/2,0), (1/2,0), and (0,1)^) with 0 < < V^f2, and —)■ 0 as i —)■ oo. The 
generic conhguration specified by (Al)-(A6) is depicted in Figure 4.2. For any sequence of polygons satisfying 
our assumptions, we now show that the interpolation error cannot be bounded by the iL^-semi-norm for a 
particular function. 

Lemma 4.3. For a sequence of polygon-vertex-edge tuples (Ti, ej)“ satisfying (Al)-(A6), 

\u — Ip.u\jji,p.^ 

lim -j—j- = oo, 


for the function u(x, y) = x'^. 

Proof. First note that |u|p 2 (p.) is trivially bounded by a constant: 

(4.4) |m|p-2(p.) = j 2^ = 4 |Pi| < tt, 

where the final inequality follows from the isodiametric inequality [37, p. 69]. 

It remains to be shown that \u — Ip.u\ pi grows without bound as i —>■ oo. We will show that the 
partial derivative with respect to y blows up near the origin. Let hi(x) denote the height of Pi at any value 
of X along edge e^, as shown in Figure 4.2. The convexity of Pi, along with assumptions (A4) and (A5), 
ensures that the region bounded by —c„ < x < 0 and 0 < y < /ii(x) is a subset of Pi. We will integrate over 
the even smaller domain, where —c^/2 < x < 0; note that assumption (Al) ensures c„ < 1. For u = x^. 
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Fig. 4.2: Generic configuration of a convex polygon of diameter 1 assumed for Lemma 4.3. For u{x, y) := x'^, 
it is shown that integrating the partial derivative with respect to y of the interpolant of u over the shaded 
region grows without bound as d{vi, e^) —)■ 0. 


= 0 and thus 


(4.5) \u- IpMm{p^)> 1^ = L - 


dy 


f^hi (fc) 


f-cl/2 Jo 


—Ip.u(x,y)] dydx. 


The inner integral is the iL^-semi-norm of a harmonic function, so Dirichlet’s principle provides a lower 
bound. Using Corollary B.2 from Appendix B, we have 


(4.6) 


fhiix) ^ Q 

—Ip^u{x,y)] dy> 


{Ip^u{x,hi{x)) - lp-u{x,0)) 

hi{x) 


We next show that lp.u{x,0) — Ip-u{x, hi{x)) > c^/2. First note that by assumptions (A3) and (A5), edge 
Ci contains the interval [—Cy,Cy] lying on the x-axis. Since lp.u{x,0) linearly interpolates the values of u at 
the endpoints of edge and u = x^ hy hypothesis, we have 


(4.7) 


lp-u{x,0) > cl- 



Fig. 4.3: On [—1,0], any piecewise linear interpolant of x"^ lies above x"^ and below |a;|. 
provides an estimate for Ip.u{x^ hi(x)) in the notation of in the proof of Lemma 4.3. 


This simple fact 


Next, note that for x G [—c^, 0], the one-variable function Ip.u{x, hi{x)) reads off the data along a portion of 
the boundary of Pi. In particular, it is a piecewise linear approximation of x'^ that equals x"^ at the x values 
corresponding to the cc-coordinate values of the vertices of Pi. We have —1 < —Cy by assumption (Al) and 
lp-u{0, h{0)) = 0 since u = x^. Thus < Ip.u{x, hi{x)) < |a:| as shown in Figure 4.3. We conclude that for 
< a: < 0, 


(4.8) 


Ip.u{x,hi{x)) < c^/2. 







By (4.7) and (4.8), we have lp-u{x,0) — Ip.u{x,hi{x)) > c^/2, as claimed. By assumption (A4), hi(x) < 
dist(vj,ei) for x S [—cJ/2,0]. Thus, the right-hand side of (4.6) is bounded below by (c„/2)^/dist(vi, e^). 
Putting this together with (4.5), we have that 


2 ( 




dx = 


_c2/2 dist(vj,ei) 8dist(vi,ej) 


Finally, using (4.4) we have that 


1 '“ ^Pi'^\H^(Pi) 


> 


lH2(p,) 


87r dist(vi, e^) ’ 


which grows without bound as i —>■ oo by assumption (A6). □ 

4.3. Example: A Prototypical Family of Non-convex Polygons. Before attempting to expand 
the result of Lemma 4.3 to the case of non-convex polygons, we will consider a simple example containing 
representative non-convex geometric configuration. In the example, we begin by attempting to bound the 
interpolation error despite the fact that most of the assumptions of Section 4.2 hold. The analysis connects 
success of the interpolation estimate to the existence of a certain discontinuous function on the boundary 
of the domain in a particular Sobolev space. For the convex polygons in question, the particular Sobolev 
space on the domain boundary) does not contain discontinuous functions and thus the estimate fails. 

When considering three-dimensional polyhedra, however, the Sobolev embedding does allow the suitable 
discontinuous functions and thus this example provides a template for the main result in Section 5.3. 




Fig. 4.4: (a) a family of non-convex polygons satisfying assumptions (A1)-(A3) and (A5)-(A6), which is 
analyzed in Section 4.3. (b) The discussion makes use of the domain which is a subset of as shown. 


Let Pe be the five-sided non-convex polygon with vertices at (—1,0), (—1,1), (1, 0), (1,1), and (0, e), with 
0 < e < 1, as depicted in Figure 4.4a, and consider interpolating the function u{x,y) = on P^. We will 
use the notation Pq to refer to the limiting domain of the union of two triangles. Observe that the sequence 
associated to {Pi/i} for i = 2,3,4,..., when scaled by a factor of l/\/5, satisfies assumptions (A1)-(A3) and 
(A5)-(A6). The scaling by the fixed value l/-\/5 does not affect the subsequent argument, so we focus on 
estimating \u — / \u\p 2 (^p^) as e —>■ 0 without explicitly scaling the domain to diameter one. Note 

also that the constrained Delaunay triangulation of P^ contains a triangle with vertices (—1,0), (1,0) and 
(0, e) for which the circumradius grows as e —> 0, so Theorem 4.1 cannot be employed. 

Observe that |w|^ 2 (p ^ is bounded below uniformly and \u\pi^p j is bounded above uniformly: 

(4.9) > |u|p2(pg) = f 2^ = 4; 

Pq 

(4.10) \'^\h^(p,) — \'^\h^(p-i) ~ ( (2^^) = 8/3. 

J Pi 
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By the triangle inequality 
(4.11) 


thus the remaining term to estimate is \Ip^u\^xf^p y Standard approaches toward this goal (e.g. the Sobolev 
embedding theorem, Morrey’s inequality) cannot be established uniformly in a trivial way due to the non¬ 
locality of the boundary, where the domain becomes arbitrarily narrow near the origin. To skirt this difficulty, 
we work with the subdomain := {(x,y) € P^l x > 0}, as shown in Figure 4.4b. 

Since u{x,y) = x'^, the boundary condition gg : dP^ —>■ M is identically 1 along the horizontal and 
vertical edges of Pg and decays linearly from 1 at (±1,1) to zero at (0, e) along the two diagonal edges. 
By Dirichlet’s principle (see Theorem B.l from Appendix B), Ip^u has minimal iJ^-semi-norm among all 
functions satisfying this boundary condition. By the symmetry of the domain and the boundary condition, 
Dirichlet’s principle also holds over the subdomain Qg, i.e., 

( 4 - 12 ) = 2 \IPtU\pif^Qy < 2 , 

for any v G that satisfies the same boundary condition as Ip^u on all the boundary segments shared 

between Pg and Qg. In particular, (4.12) applies in the case v := Tr“^(jig, where Tr“^ denotes the right- 
continuous inverse of the trace operator for the domain Qg. Note that the trace operator and its inverse 
depend on the domain, i.e. on e, but we omit this explicit dependence from our notation as dependence on 
e is already indicated by the boundary data ^g. We can now estimate 

(4.13) l^f’e'*^lpi(Pg) — 2 |Tr 9e\pi(^Q^^ ^ 1^4.13 \ \9e\\H^/^(dQe) ' 

For Lipschitz domains, the trace operator has a well-defined continuous, bounded inverse, meaning that for 
any fixed e, the constant C above exists; see Theorem A.2 in Appendix A for more details. 

Before analyzing ||ffe||pi/ 2 (aQ^), we must argue that the constant 04,13 is independent of e. The typical 
technique in the analysis of Sobolev-type properties on Lipschitz domains, e.g. extension or trace theorems, 
involves covering the domain by a finite number of overlapping patches, constructing a partition of unity 
supported by the patches, and then on each patch intersecting and “flattening” the boundary via a Lipschitz 
map. The bound on the inverse trace operator depends on the domain only by the parameters of this 
decomposition: the number of mutually overlapping patches, the behavior of the partition of unity, and the 
Lipschitz constant of the maps for flattening the boundary. For Qg, this construction can be performed 
with a set of patches and partition of unity that is independent of e. Further, the Lipschitz constant of the 
flattening maps can also be bounded independent of e as shown in Lemma A.3 in Appendix A. Thus the 
constant < 74.13 can be taken independent of e. 



Fig. 4.5: Contour plots of Ip^u for (left to right) e = 0.5, e = 0.1, and e = 0.02. Large gradients are 
concentrated near a single point, but, in two dimensions, this is not not enough to cause the interpolation 
error to be unbounded in the limit. 


We examine the limiting behavior of the estimate by considering the limit of go, i.e. the limit of the 
boundary condition. Here, go is discontinuous on the one-dimensional boundary of the limiting two-triangle 
domain Qo- In one dimension, the space contains only continuous function, as guaranteed by the 

Sobolev embedding theorem or, more specifically, by Morrey’s inequality (stated in Theorem A.l in Ap¬ 
pendix A). Thus ge does not converge in as e —>■ 0 and hence || 5 e||/ 4 i/ 2 (aQ^) cannot be bounded 

uniformly. 

Remark 4.4. The Sobolev embedding theorem ensures the continuity of functions in when 

mp > d where d is the dimension of the space. The setting relevant to the example above (d = 1, m = 1/2, 
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p = 2) is the equality case of this inequality. The case of polyhedra in three dimensions involves combinations 
that do not satisfy the inequality (d = 2 , to = 1 /2, p = 2) and the resulting conclusions will be different. 
These cases are discussed in Section 5.3. Similarly, analyzing harmonic coordinates in spaces for p > 2 

can also lead to error bounds, even for polygons; in those cases, however, there are other hurdles, such as 
the inability to directly appeal to Dirichlet’s principle. 

4.4. Sharpness: Non-convex Polygons. To analyze the sharpness of estimate (4.1) for families of 
non-convex polygons, we begin with the same assumptions that preceded Lemma 4.3, except for (A4), which 

is essentially tied to the local convexity of the polygon at the vertex in question v^. This will be replaced 

by some other restrictions to ensure that non-adjacent entities of the polygon do not interfere with the 
interaction of vertex and non-incident edge e^. 

Specifically, we restrict ourselves to families of polygons satisfying (A1)-(A3), (A5)-(A6), and the fol¬ 
lowing: 

(A7) dist(vi,Va) > Cy for all vertices of Pi where ^ v^; 

(A8) dist(vi,ea) > c„ for all edges of Pi other than and the two edges incident to v^; 

(A9) the segment between and the origin is contained in Pi. 

In the above assumptions, c„ is the same constant as in (A5), although this holds without loss of 
generality since it can be assumed to be the minimum constant among (A5), (A7), and (A8). For polygons 
with at least four sides, (A8) is strictly stronger than (A7) since any vertex near also will be connected to 
an edge that is not incident to v^. Still, we include (A7) for clarity. The assumption (A9) serves to exclude 
situations in which approaches from the outside of the Pi which generally does not cause the estimate 
to blow up; see Figure 4.6. 





Fig. 4.6: Several non-convex polygons that violate various restrictions in Lemma 4.5: (A) a polygon with 
two nearby vertices, violating (A7); (B) a polygon with a vertex near multiple non-incident edges, violating 
(A8); (C) a polygon with a vertex approaching a non-incident edge from the outside of the polygon,violating 
(A9). 


Lemma 4.5. For a sequence of polygon-vertex-edge tuples (Pi, I’i, satisfying (A1)-(A3) and (A5)- 

(A9), 


\u 

lim -^^- = (X), 

2—)-00 


for the function u{x, y) = x^. 

Proof. First, we identify some some implications of the assumptions on the domain and, without loss of 
generality, define a generic case for analysis, which is depicted in Figure 4.7. Specifically, let di^i and di ^2 be 
the two edges of Pi incident to vertex Vi. Without loss of generality, we assume that di^i lies in the second 
quadrant, i.e. a: < 0, since reflecting the domain across the y-axis does not change the subsequent analysis. 
Moreover, if di ^2 is also in the second quadrant, we assume that di^i is “below” di^ 2 , or, more formally, 
a counter-clockwise sector from i to di ,2 is interior to the polygon. Finally, we restrict our analysis to 
the case where the line through di^i has negative slope: if the line has positive slope, then we have the 
configuration as in the proof of Lemma 4.3 and can apply the same argument. 
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-(c.)2/V2 

Fig. 4.7: Generic configuration of a polygon assumed for Lemma 4.5, showing the standard and rotated 
coordinate systems {x^y) and (x,y), respectively. For u(x,y) := it is shown that integrating the partial 
derivative with respect to y of the interpolant of u over the shaded region grows without bound as d{vi, Ci) —)■ 

0 . 


As in Lemma 4.3, we will identify a suitable integration subdomain and show that the interpolation 
estimate must blow up simply looking at the subdomain. Estimate (4.4) still holds, providing a constant 
upper bound on the denominator: < tt. For the numerator, first note that 




The term can also be bounded above by a constant, similar to |n|p 2 (p.), and thus we only need to 

show that |-fPiM|pi(p.) grows without bound as i ^ oo. 

Establishing a lower bound on \IPiU\H^{Pi) is similar to (4.5), but involves analyzing a different integration 
domain. To describe this domain, consider an alternative coordinate system, with axes denoted x and y, 
centered at and rotated 45° clockwise: in other words, the x axis lies on a line of slope —1 and the y 
axis lies on a line of slope +1, as shown in Figure 4.7. We restrict attention to the region between edges 
and di^i in the range x G (—c„/2, 0), which is the shaded region in Figure 4.7. Denote this region by Dp 
For sufficiently large i, assumption (A6) in conjunction with (A7)-(A9) ensure that Di C Pi] from here on, 
we always take i large enough that this containment holds. Further, the {x, y) coordinate system has been 
constructed so that in our generic configuration, with dj i downward sloping in the second quadrant, the 
distance between di^i and is small in the y direction. Estimating |7p;U|pi^p.^ over the domain Di gives: 


(4.14) 


l-^Pi'^lpi(Pi) ^ 


'Di 



> 


f / 

J-cl/2 JLiix) 



IPiU{x,y) 


2 


dy di:, 


where Di represents the integration subdomain associated with polygon Pi and Li {x) denotes the appropriate 
slice of the integration domain for a given x value. To estimate the inner integral, we observe that Dirichlet’s 
principle can be applied as was done in (4.6) and thus. 



d {Ip^u{x,di^i{x)) - Ip.u{x,ei{x))f 

ip.P.ui..y)) iy> -- 


where {x,di^i{x)) and (x,ei(x)) denote the end points of Li(x) on di^i and ei(x), respectively, and \Li{x)\ is 
the length of segment Li{x). 

As in Lemma 4.3, u along a is at least cj, and u along di^i is at most cl/2. Thus (4.14) becomes. 


(4.15) 


I^Pi'^lpi(Pi) ^ 


-cl/2 4|T*(x)| 


dx. 


The integration set Li{x) is contained in the second quadrant and a line of slope +1, and thus we can 
estimate 


\Li{x)\ < 2 {dist{-vCi)/V2 - x^ 
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noting that x < 0 in our area of interest. Thus (4.15) becomes 


\^Pi'^\H^{Pi) — 


Tc 2/2 8 (dist(vj,ei)/v^-x) 8 


dx = y (in (^dist(vi,ei)/y 2 + c^/ 2 ^ - In (^dist(vi, ei)/y 2 ^^ 


Taking the limit as i —>■ oo, the first term on the right approaches a constant while (A 6 ) ensures that the 
second term grows without bound, thus completing the result. □ 

5. Harmonic Coordinates on Polyhedra. As is typical in computational geometry, going from two 
to three dimensions introduces additional subtleties and challenges. On tetrahedra in as on any n-simplex 
in R”, barycentric functions are uniquely defined and linear, yielding a similar analysis of their interpolation 
properties. On a general polyhedron P, however, just defining generalized barycentric coordinates on dP is 
not unique or even trivial, since for non-triangular boundary facets, a two-dimensional generalized barycentric 
coordinate must be employed. Once values on dP are hxed by some means, they are used as boundary 
conditions to define the harmonic coordinates by the same PDE (3.1) from the two-dimensional construction. 

In Section 5.1, we restrict our analysis to polyhedra with triangular facets, where the GBCs on the 
boundary are linear functions. This is the most common setting for existing work on 3D GBCs [40, 48, 46, 47] 
although some constructions have been developed and analyzed more generally [89, 90, 44, 39]. While the 
result should still hold for a broader class of polyhedra with non-triangular faces, the restricted setting 
helps to avoid complexity associated with the two-dimensional GBCs on the boundary facets, for which 
shape-quality and error estimates are already signihcantly more intricate than in the triangular setting. 

In Section 5.2, we consider arbitrary convex polyhedra and demonstrate that the error estimate will fail 
in any family of polyhedra where vertices are allowed to be arbitrarily close to the interior of non-adjacent 
facets. Finally, in Section 5.3, we show that non-convex polyhedra do not have this property; vertices can 
approach opposite faces in certain ways such that the error does not blow up. Comparisons to the 2D 
analogues of each result are made in each section. 

5.1. Convex, Bounded Aspect Ratio Polyhedra. First we dehne the set of polyhedra for which 
error estimates analogous to the result in [41] can be established. For a convex polyhedron P, an insphere of 
P is a sphere inscribed in P of maximum radius p(P). Given P, the inradius p{P) is uniquely defined and 
we fix some particular insphere center c € P arbitrarily if there is not a unique insphere. The inradius of a 
two-dimensional facet P of P is defined as in the polygonal case and denoted p{F). The aspect ratio, also 
called the chunkiness parameter, is denoted by 7 and is defined as the ratio of the diameter to the inradius 
of P, i.e. 


diam(P) 

P{P) ■ 

Given a bound 7 * > 0 on the aspect ratio of both the polyhedron and each of its faces, let V be the set 
of polyhedra P with the following properties: 

• P is convex; 

• all facets of P are triangles; 

• l{P) < 1*\ 

• p{F) > diam(P)/ 7 * for each facet F of P. 

These restrictions immediately impose a bound on the number of faces and vertices of the polygon as 
stated in the following lemma. 

Lemma 5.1. There exists > 0 depending only upon 7 * such that all polygons ofV contain fewer than 
n* triangles and vertices. 

Proof. Without loss of generality (up to affine scaling), we assume that P has diameter 1. Since the 
filled sphere has maximum surface area over all convex sets of equal diameter [ 11 ], the surface area of any 
P G P is at most tt. The inradius of each face of P is at least -P, so estimating very conservatively, there 
are at most TryJ triangles forming the boundary of P. Fuler’s formula for the boundary triangulation of the 
polygon gives the following relationship between the number of vertices and the number of triangles np. 
ny = nt/2+ 2 < nt, since nt is at least four. Thus both the number of triangles and vertices on the boundary 
of the polygon are bounded by n* = TryJ. □ 
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On the class of bounded aspect ratio polyhedra with triangular facets of bounded aspect ratio, the 
standard interpolation error estimate holds. 

Theorem 5.2. There exists a constant C > 0 depending only upon 7* such that for any polyhedron 
P G V and all functions u G H^{P), 

(5.1) |u - < C diam(P) |m|^2(p) • 


Proof. The proof consists of two main ingredients. First we observe that the majority of the proof 
of Theorem 4.1 generalizes to this setting, leading to the conclusion that the interpolation error in the 
harmonic interpolant is no worse than that of the piecewise linear interpolation using any tetrahedralization 
that respects the triangle-faceted boundary. Second, we show that a bounded aspect ratio tetrahedralization 
of any polygon in V can formed using a single additional interior vertex at the incenter. This is the same 
argument used in for the polygonal case in [41, Theorem 2]. 

As usual, we restrict our analysis to polyhedra with unit diameter, without loss of generality. Following 
the proof of Theorem 4.1, applying the Bramble-Hilbert lemma and Dirichlet’s principle, we first estimate 
the interpolation error by 

\u — Ipu\pi^p^ \u — Piilpi(p) + \Pu ~ ^p''^\h^(p) 

^ \u — Pu\pi(^p'j + \It {Pu — R)lgifT1 

tgt 

(5.2) < 2 ju — Pujpi^p^ + ju — Ipulpifp) . 

T&r 

As before, is the linear polynomial from the Bramble-Hilbert lemma, but T can be any tetrahedralization 
of P, possibly including additional interior vertices. We require that each face of P must be the face of one 
of the tetrahedra in T, i.e. there can be no mesh refinement along the boundary. 

It remains to establish an estimate on ju — As the standard interpolation error for the linear 

interpolant on tetrahedra, this term can be estimated depending only upon the aspect ratio of T. Select T to 
be the tetrahedralization formed by adding an additional vertex at an incenter of P and forming tetrahedra 
by connecting this incenter to each of the faces. For each tetrahedron T G P, a, rigid body transformation 
will place the face of T forming part of the boundary of P in the ccy-plane. 

By convexity, the aspect ratio bound on V and the definition of the incenter, the distance from the 
incenter to the xy-plane (i.e., /i* in Lemma E.l) is at least I/ 7 *. Since V requires the inradius of each face 
to be at least I/ 7 *, Lemma E.l asserts the existence of a quality bound yp depending only upon the original 
aspect ratio bound 7 *. 

The standard interpolation error estimate on tetrahedra asserts the existence of a constant Cp depending 
only upon the aspect ratio of T such that ju — Ipujpi^p) < C'pdiam(r) |m|p 2 (p). Applying this to (5.2) and 
using the Bramble-Hilbert estimate completes the result: 

\u — Ipu\hi(^p^ < 2 Cp |u|p2(p) = {2 + Cp) |w|p2(p) ■ 

Ter 


□ 


5.2. Convex Polyhedra Without Bounded Aspect Ratio. When translating the aspect-ratio in¬ 
dependent error estimates on convex polygons in Sections 4.1 and 4.2 to three dimensions, much of the 
analysis is identical albeit without the connections to constrained Delaunay tetrahedralization. The inequal¬ 
ities (5.2) hold under a very limited set of assumptions, meaning that, in general, harmonic coordinates 
always perform at least as well as using a piecewise linear interpolant over a tetrahedralization of P. Con¬ 
structing an analog to Theorem 4.1 is not so straightforward, since interpolation error on tetrahedra is not 
bounded by circumradius, due to a class of tetrahedra called slivers. Slivers are nearly coplanar tetrahedra 
with no short edges, all four vertices near a common circle, and circumradii proportional to the edge lengths. 
These tetrahedra have poor interpolation properties, since they allow for disjoint edges that pass arbitrarily 
close to each other; see Figure 5.1. Sliver tetrahedra are not eliminated by the Delaunay construction and 
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Fig. 5.1: A sliver tetrahedron depicted from three different angles with its (modest sized) circumsphere. 


are a well known detriment to 3D Delaunay meshes; a substantial body of mesh generation literature is 
devoted to sliver removal, e.g. [23, 33, 32, 59]. 

Despite a great deal of interpolation error estimates for tetrahedra without a bound on aspect ratio 
(e.g. [45, 52, 76]) the general error metric, called “coplanarity” in [70], does not correspond to a natural 
geometric construction. While Theorem 4.1 can be established in terms of a coplanarity measure, there is 
no apparent tool for generating tetrahedralizations of a given polyhedron with minimal coplanarity, even if 
it is assumed that all faces are triangles of high geometric quality. Worse still, some polyhedra cannot be 
tetrahedralized without the introduction of additional vertices [56, 75] including some convex polyhedra [68]. 
Accordingly, the constrained Delaunay tetrahedralization is not well-defined in general and typical approaches 
to ensure its existence involve adding additional vertices on the boundary [77] , which would interfere with the 
construction leading to estimate (5.2). While tetrahedralization does not provide a sharp characterization 
of interpolation error, many of the general principles still hold. In particular, looking at the sharpness result 
in Section 4.2, we now see that interpolation error is still guaranteed to grow whenever a vertex is very near 
the interior of a non-adjacent face. 

Building the analog of Lemma 4.3, consider a sequence (Pi, Vi, where Pi is a convex polyhedron, 

Vi is a vertex of Pi and /i is a face of Pi, under the following assumptions: 

(Bl) diam(Pi) = 1; 

(B2) fi lies in the xy-plane; 

(B3) Vi lies on the positive z-axis; 

(B4) There exists a constant c„ such that dist(o, 9/i) > where dfi is the boundary of ft and o is the 
origin; 

(B5) dist(vi, /i) —>• 0 as i —>■ oo. 



Fig. 5.2: Example convex polyhedron containing a vertex (black dot) which is near the interior of a non- 
adjacent face. The integration domain where the gradient is known to be large is depicted with the dark 
gray half circle. 
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In this setting, the following lemma asserts that the standard error estimate cannot be established. 
Lemma 5.3. For a sequence of polyhedron-vertex-face tuples (^i, satisfying (B1)-(B5), 


,. 1'^ 

lim - 1 — I -— - = oo, 


for the function u{x, y, z) = x^ -\- y'^. 

Proof. The proof follows the same steps as the 2D case, Lemma 4.3. The isodiametric inequality again 
implies a uniform upper bound on the denominator, reducing the desired estimate to showing that the 
numerator grows without bound. 

Again, the integral in the numerator is restricted to the key portion near the origin, where large gradients 
occur in the z-direction: 




pO |^hi{x,y) / Q \ 2 


Above, hi{x,y) denotes the height of the boundary of the polygon above the cry-plane (which is well-defined 
near the origin above fi). Here Di is the integration domain in the acy—plane: half of the ball centered 
at (0,0) with radius cl/2, where the plane supporting Pi at (0,0,/ii(0, 0)) is below hi{0,0); see Figure 5.2. 
The inner integral is estimated using the piecewise linear nature of the interpolant on the boundary and we 
conclude that, in the integration domain, 

Ip^u{x,y,hi{x,y)) < c^/V2. 


Then after applying Corollary B.2 to estimate the inner integral, we can estimate the remaining terms: 


1^ — 


iDi 


{Ip.u{x, y,hi{x, y)) - Ip,u{x, y, 0)) 
hi{x) 


<lA{x,y) 


> 


r [( 1 - 1 /V 2 ) eg 

Id, dist(vi,/i) 

TTC^ ^ [(1 - 1/72) eg]' 
8 dist(vj,/i) 


■<lA{x,y) 


As dist(vi, fi) is the only remaining non-constant term, (B5) ensures that the estimate grows without bound. 

□ 

5.3. Non-convex Polyhedra. Non-convexity presents additional challenges when moving to three 
dimensions. Remark 4.4 suggested that vertices of non-convex polyhedra can approach non-adjacent faces 
without causing the interpolation error to blow up. The lemma below formalizes the three dimensional 
analog of the example in Section 4.3 demonstrating the a sharp characterization of polygons that admit 
bounded interpolation error estimates must include some shapes that are fundamentally different than those 
allowed in two dimensions. 

Lemma 5.4. There exists a sequence of polyhedron-vertex-face tuples [ Pi , Vi, fi)/Z^, satisfying (B1)-(B5) 
such that 


lim 
2—>-00 


W\h^{P,) 


= C < oo, 


for the function u{x, y, z) = x'^ -\- y^ and some real number C. 

Most of the details are are a direct generalization of the construction given in Section 4.3 to three 
dimensions. Here we outline the key points and differences from two dimensions. The construction involves 
polyhedra shown in Figure 5.3a with 13 vertices: (1,1,0), (—1,1,0), (1, —1,0), (-1,-1, 0), (1,1,1), (—1,1,1), 
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(b) 


Fig. 5.3: (a) a non-convex three-dimensional analog to the domains in Figure 4.4a used in Lemma 5.4. The 
family of polyhedra involves allowing the vertex at the interior of the convex hull to approach the opposite 
face, (b) The convex quarter domain needed for analysis. 


(1,—1,1), (—1,—1,1), (1,0,1), (0,1,1), (—1,0,1), (0, —1,1) and (0,0, e). The vertex (0,0, e) is allowed to 
approach the non-adjacent face in the xy-plane. 

If all faces of these polyhedra are triangular, then (0,0, e) is approaching a point on the boundary of a 
triangle and thus (B4) is not satisfied. This can be corrected by adding vertices to the face of the polygon in 
the xy-plane so that (0, 0) lies in the interior of a triangle. We ignore this minor detail, allowing polyhedra 
with non-triangular (in this case square) faces. 

In the analysis, one fourth of the original domain is considered, namely, the portion of domain with 
a: > 0 and y > 0 shown in Figure 5.3b. Symmetry of u across the xz- and j/z-planes ensures that analysis 
on this subdomain is sufficient: any integrals on the subdomain will be one fourth of the integral over the 
entire domain. As in the family {Qe} in Section 4.3, the resulting domains are convex and have a similar 
“uniform Lipschitzness.” 

Rather than elaborate upon the details, we emphasize the key properties of the limiting case that are 
driving the result. The limiting boundary conditions have a discontinuity at a single point (the origin), 
but, in two dimensions, discontinuous functions are admitted in including these specific boundary 

conditions. Thus, the limiting interpolant Ip„u is in and is in essence providing a uniform bound on the 
iL^-norm of interpolants for the cases with e sufficiently near 0. 

The conclusion to be drawn from Lemma 5.4 is as follows: the class of non-convex polyhedra that admit 
bounded interpolation errors is fundamentally different than all of the other geometry classes considered in 
this work. For simplices (in any dimension), polygons (convex or non-convex) and convex polyhedra, a se¬ 
quences of domains involving a vertex approaching an opposite edge/face causes an unbounded interpolation 
error while Lemma 5.4 demonstrates that this is not the case for non-convex polyhedra. 
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Appendix A. Sobolev Spaces. 

For multi-index a = ( 01 , 02 ) and point x = {x,y), define x“ := o! := 01 O 2 , |o| := oi -I- 02 , and 

:= ujdx°‘^dy°‘'^. The Sobolev semi-norms and norms over functions defined on an open set are 

defined by 

Hwrr^.rin)-= f XI l^““Wrdx and := ^ 

\a\=m 0<k<m 

In the case p = 2, these norms define Hilbert spaces and a denoted := The case m = 0 

reduces to space of integrable functions: the common spaces. 

We now state a particular Sobolev embedding property that is most relevant to our analysis; this result 
is phrased in much more generality in the literature on Sobolev spaces, e.g., [3, 57]. 

Theorem A.l. (Morrey’s inequality) If ill € K” is a Lipschitz domain, and m > n/p, then all functions 
in lT’”’^(r2) are continuous. 

This theorem is valid for non-integral m, i.e. for fractional Sobolev spaces [29]. In our context, these 
fractional Sobolev spaces arise in the definition of trace spaces: the space of functions corresponding to 
the boundary values of the a Sobolev function. The trace theorem below asserts that such functions are 
well-dehned in the simplest setting. 

Theorem A. 2. (Trace Theorem) [63, 26, 30] If Tt £ K" is a Lipschitz domain, then the trace operator, 

Tr : H^{n) H^/^{dVl) 

is a well-defined, hounded, linear operator with a bounded right-inverse, i.e. there exists a constant C > 0 
such that for all u € and for all g S H^/'^(dTi), 

||Tr 'u||//i/ 2 (af 2 ) < ^ 1iRi(n) 

^5|liyi(n) — ^ll5llHi/2(an) ■ 


The typical theory of Sobolev spaces on Lipschitz domains involves covering the domain by patches and 
flattening the boundary on each patch to transform the domain locally into a half-space. The following lemma 
is a technical detail used in Sections 4.3 and 5.3 when arguing that certain domain-dependent constants are 
in fact domain independent for a specific family of domains. 

Lemma A. 3. For any e G (0,1/4), the domain (defined in Figure 4-4) can be locally flattened using 
decomposition with a uniformly bounded Lipschitz constant that does not depend on e. 

Proof. Consider flattening the domain in three different patches shown in Figure A.l. The origin of 
each of the local coordinates systems is located at one of the vertices of that does not depend on e. We 
will show that the Lipschitz transformation for the patch centered at the origin is uniformly bounded. The 
patch centered at (1,1) is similar and the patch centered at (1, 0) is simpler because locally the domain 
is independent of e. 

The local coordinate system is defined by the following (orthogonal) transformation: 


1 2 

X = —X - 


75 

2 

V = —;=^ 

75 


75^ 


In this local coordinate system and for a ball of at radius at least 1/2, the location of boundary of can 
be written as a function b of x: 


b{x) 


2 x 
— Ix 

3-e y. I VEe 
. 2e-l‘^ 2e-l 


if i > 0; 

if - ^ ^ ^ 7 

ifi<-^e. 
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Fig. A.l: Local axes used for flattening the boundary of the domain (see Figure 4.4 for context). These 
local axes and patches can be used for all e < 1/4. 


For e < 1/4, we have < 11/2, so the derivative of b is uniformly bounded with respect to e. □ 

Appendix B. Dirichlet’s Principle. 

Dirichlet’s principle asserts that harmonic functions have minimal iL^-norm among all functions with the 
requisite boundary conditions. A precise statement and proof in a general setting can be found in Evans [36, 
Chapter 2, Theorem 17]. For reference, we state Dirichlet’s principle in the notation of this paper. 
Theorem B.l (Dirichlet’s principle). Assume Ipu € C^(P) solves (3.3) and let 

A\= {ic G C‘^{P) : w = gu on dP}. 


Then 


(B.l) 


\Ipu\hi(p) = mn|w|^i(p) . 


Conversely, if Ipu G A satisfies (B.l), then Ipu solves (3.3). 

In ID, the harmonic interpolant is the line fitting the values at the ends of an interval. In this setting, 
Dirichlet’s principle can be given as a more concrete inequality. 

Corollary B.2. Let Ca,Cb G K and let /* : [a,b] M. be the linear function with fstfa) = Ca and 
f*{b) = Cb. Then for all f G H^{[a, 6 ]) such that f{a) = Ca and f{b) = Cb, 




da; > 




da; 


jcb - CaY 

b — a 


Appendix C. Bramble-Hilbert Lemma. 

The Bramble-Hilbert lemma [15] is the key estimate approximating Sobolev functions by polynomials in 
the appropriate context for the finite element method. While the basic finite element analysis only requires 
the theorem to hold on a specific reference element, the theorem can be established uniformly over convex 
domains [85, 28] which is essential to the analysis over general polygonal/polyhedral domains. 

Theorem C.l. There exists a constant Cc.i depending only upon n such that for any convex polytope 
of unit diameter P C K” and any u G H^{P), there exists an affine function such that 

(C.l) \u — Pu \ pl (^ p ^ C : Cc . l \ u \ jj 2 (^py 
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Fig. D.l: A set of points and constrained segment for constrained Delaunay triangulation. The triangle shown 
is not part of the Delaunay triangulation of the points but is in the constrained Delaunay triangulation, since 
all vertices inside its circumcircle are not visible through the constrained segment. 



Fig. D.2: In a Delaunay triangulation, the triangle adjacent to an obtuse triangle along the (unconstrained) 
longest edge always has a larger circumradius. 


Appendix D. Constrained Delaunay Triangulation. 

Given a set of points V = {v^}, the Delaunay triangulation consists of the set of triangles T = {Tj} 
formed from vertices of V satisfying the empty circumball property: the circumcircle of Tj contains no points 
of V in its interior. This uniquely defines a triangulation when the points of V lie in general position, i.e., 
no four points share a single circle, but can be defined for any point sets using some arbitrary tie-breaking 
rules. 

Given a set of non-intersecting segments S = {<5'^} where each segment has end points in V, the con¬ 
strained Delaunay triangulation [24, 77] is defined with a modified empty circumball criteria to create a 
triangulation with some similar properties while ensuring that the segments of S appear in the triangula¬ 
tion. Specifically, the empty circumball property is relaxed by allowing points in the triangle’s circumball if 
they are not “visible” from the triangle. A point is visible to a triangle if every line between that point and 
any point in the triangle does not cross a required segment; see Figure D.l. 

Lemma D.l. Let Ti and be adjacent triangles in a constrained Delaunay triangulation along a non- 
constrained edge E. If Ti is obtuse and adjacent to T 2 along edge E then the circumradius ofTi is smaller 
than the circumradius 0 /T 2 . 

Proof. Let v denote the vertex of Ti opposite edge E and let Ci and C 2 denote the circumcenters of Ti 
and T 2 , respectively; see Figure D.2. Since Ti is obtuse, Ci lies outside of Ti. Since the circumcenters are 
equidistant from the endpoints of E, Ci and C2 lie on the line perpendicular to E. If C2 lies on the side of 
this line nearer to Ti, then the circumcircle of T 2 contains v and thus is not in the constrained Delaunay 
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triangulation. So C 2 lies on the side of the line away from Ti and thus the circumradius of is larger than 
that of Ti. □ 

Appendix E. Tetrahedron Quality. 

Fixing constants r* > 0 and > 0, let T 2 denote the set of triangles in the xy-plane with inradius 
larger than r* and diameter at most 1. Also define a set of points 

V = {{x,y,z): + j/^l < 1, 0 << z < 1 } . 

Let Ts denote the set of all possible tetrahedra with one face a triangle from 72 and the fourth vertex from 

V. 

Lemma E.l. There exists a constant 7 ^ < 00 depending only upon r* and h* such that the aspect ratio 
of every tetrahedron in Ta is smaller than ^t- A sketch of the proof of this lemma is as follows. Since 
r* > 0 and /i* > 0, the aspect ratio of every tetrahedron in V is strictly larger than 0. Consider a sequence 
of tetrahedra {Ti) such that limit of the aspect ratio approaches the infimum over V. Via compactness, we 
can extract a subsequence of this sequence such that the vertices converge. Thus by our original statement 
(all elements of V have strictly positive aspect ratio), the lower bound on the aspect ratio of all tetrahedra 
in V is strictly positive. 
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